{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "81951160-fcd9-43ab-b9ee-b990b6c3606c",
   "metadata": {},
   "source": [
    "# Get metrics: heatwave, polluted day and wildfire"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b0a650fa-143d-4cb5-8ab8-0a3be1142887",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 800x600 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import helper_400\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "helper_400.set_sns_style()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "dac02566-bd55-4971-ade0-6f00e0c38e4e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = pd.read_parquet(\"outputs/final_merge_5_western_us.parquet\", engine=\"pyarrow\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "629d3e76",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>wfday</th>\n",
       "      <th>tmax</th>\n",
       "      <th>smoke_pm</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>time</th>\n",
       "      <th>location_label</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">2006-01-01</th>\n",
       "      <th>04001942600</th>\n",
       "      <td>0.0</td>\n",
       "      <td>15.978498</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>04001942700</th>\n",
       "      <td>0.0</td>\n",
       "      <td>14.836937</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>04001944000</th>\n",
       "      <td>0.0</td>\n",
       "      <td>14.192923</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>04001944100</th>\n",
       "      <td>0.0</td>\n",
       "      <td>14.254541</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>04001944201</th>\n",
       "      <td>0.0</td>\n",
       "      <td>15.200634</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           wfday       tmax  smoke_pm\n",
       "time       location_label                            \n",
       "2006-01-01 04001942600       0.0  15.978498       0.0\n",
       "           04001942700       0.0  14.836937       0.0\n",
       "           04001944000       0.0  14.192923       0.0\n",
       "           04001944100       0.0  14.254541       0.0\n",
       "           04001944201       0.0  15.200634       0.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f11b72e3-7bcc-43ac-a823-15fa210e69b7",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time</th>\n",
       "      <th>wfday</th>\n",
       "      <th>tmax</th>\n",
       "      <th>smoke_pm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>99213732</td>\n",
       "      <td>9.921373e+07</td>\n",
       "      <td>9.921373e+07</td>\n",
       "      <td>9.794617e+07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>2013-07-02 00:00:00</td>\n",
       "      <td>1.165484e-03</td>\n",
       "      <td>2.187312e+01</td>\n",
       "      <td>5.355449e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>2006-01-01 00:00:00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>-2.898053e+01</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>2009-10-01 00:00:00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>1.574429e+01</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>2013-07-02 00:00:00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>2.205001e+01</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>2017-04-02 00:00:00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>2.851511e+01</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>2020-12-31 00:00:00</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>5.124999e+01</td>\n",
       "      <td>5.851574e+02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>NaN</td>\n",
       "      <td>3.411928e-02</td>\n",
       "      <td>9.490203e+00</td>\n",
       "      <td>4.787655e+00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                      time         wfday          tmax      smoke_pm\n",
       "count             99213732  9.921373e+07  9.921373e+07  9.794617e+07\n",
       "mean   2013-07-02 00:00:00  1.165484e-03  2.187312e+01  5.355449e-01\n",
       "min    2006-01-01 00:00:00  0.000000e+00 -2.898053e+01  0.000000e+00\n",
       "25%    2009-10-01 00:00:00  0.000000e+00  1.574429e+01  0.000000e+00\n",
       "50%    2013-07-02 00:00:00  0.000000e+00  2.205001e+01  0.000000e+00\n",
       "75%    2017-04-02 00:00:00  0.000000e+00  2.851511e+01  0.000000e+00\n",
       "max    2020-12-31 00:00:00  1.000000e+00  5.124999e+01  5.851574e+02\n",
       "std                    NaN  3.411928e-02  9.490203e+00  4.787655e+00"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.reset_index().describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a94ff5e1-a266-4ba2-b711-4479deb6f0b5",
   "metadata": {},
   "source": [
    "## Dates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f9675a5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b83d88e0-3879-4911-9220-fd11d9e113e7",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Timestamp('2006-01-01 00:00:00')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.time.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cab1124c-2c04-4789-adf6-5be79badafe9",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Timestamp('2020-12-31 00:00:00')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.time.max() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "35f473a8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = df.rename(columns={'location_label': 'GEOID'})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07e1a5cf-dec3-4547-b798-9b339715e23c",
   "metadata": {},
   "source": [
    "# Preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73862ebf-777f-4d50-9c65-480fcdd90c79",
   "metadata": {},
   "source": [
    "# Heat day - temp over 95th perc of summer months in the first 5 years as treshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "21a4c8ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time</th>\n",
       "      <th>GEOID</th>\n",
       "      <th>wfday</th>\n",
       "      <th>tmax</th>\n",
       "      <th>smoke_pm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2006-01-01</td>\n",
       "      <td>04001942600</td>\n",
       "      <td>0.0</td>\n",
       "      <td>15.978498</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2006-01-01</td>\n",
       "      <td>04001942700</td>\n",
       "      <td>0.0</td>\n",
       "      <td>14.836937</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2006-01-01</td>\n",
       "      <td>04001944000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>14.192923</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2006-01-01</td>\n",
       "      <td>04001944100</td>\n",
       "      <td>0.0</td>\n",
       "      <td>14.254541</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2006-01-01</td>\n",
       "      <td>04001944201</td>\n",
       "      <td>0.0</td>\n",
       "      <td>15.200634</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        time        GEOID  wfday       tmax  smoke_pm\n",
       "0 2006-01-01  04001942600    0.0  15.978498       0.0\n",
       "1 2006-01-01  04001942700    0.0  14.836937       0.0\n",
       "2 2006-01-01  04001944000    0.0  14.192923       0.0\n",
       "3 2006-01-01  04001944100    0.0  14.254541       0.0\n",
       "4 2006-01-01  04001944201    0.0  15.200634       0.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "dfe34de6-ff14-4079-bf46-9bc477f02899",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/lc/q1l1y0k176b0h_m5yz0s_r5w0000gn/T/ipykernel_70359/1474484182.py:11: DeprecationWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.\n",
      "  tmax95_series = df.groupby('GEOID').apply(get_tmax95_summer_first_5_years)\n"
     ]
    }
   ],
   "source": [
    "def get_tmax95_summer_first_5_years(x):\n",
    "    # Filter the data for summer months \n",
    "    summer = x[x[\"time\"].dt.month.isin([6, 7, 8, 9])]\n",
    "    \n",
    "    # Restrict to the first 5 years\n",
    "    summer_first_5_years = summer[summer[\"time\"].dt.year <= summer[\"time\"].dt.year.min() + 4]\n",
    "    \n",
    "    # Calculate the 95th percentile\n",
    "    return summer_first_5_years['tmax'].quantile(0.95)\n",
    "    \n",
    "tmax95_series = df.groupby('GEOID').apply(get_tmax95_summer_first_5_years)\n",
    "\n",
    "df['tmax95'] = df['GEOID'].map(tmax95_series)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b0e2d83a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "extreme_day_threshold = 32.22 # in C\n",
    "\n",
    "df['heatday'] = np.where(\n",
    "    (df['tmax'] >= df['tmax95']) & (df['tmax']>= extreme_day_threshold), \n",
    "    True, \n",
    "    False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "dfc72860",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAEjCAYAAABNQxmFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABpJklEQVR4nO3dd1gUx/8H8PfRVFBQEDCxxIIHoqIIdo1CFIO9axQrih17V5RorDFijTWKiAUb2HuXCLZQlCIdUToIHvXu5vcHv9vvHRz9ENDP63nyRPb2ZmdnZ/c+OzO7w2OMMRBCCCHku6ZU2RkghBBCSOWjgIAQQgghFBAQQgghhAICQgghhIACAkIIIYSAAgJCCCGEgAICQgghhIACAkIIIYSAAgJCCCGEgAICuV6/fo3FixfD0NAQPXr0wKxZs7Bw4UKMGjUKtra2uHz5coHv+Pr6olOnToiNjVVoXmbPno0//vgDAODj44OlS5fC1tZWodvIv53KduXKFWzatAn9+/eHvb09vnz5IvN5ZmYm1q1bh/Xr12PatGk4evSo3HSGDRsGQ0NDGBoaonXr1oUem9DQUDg6OmLQoEEK35eqLn+9ffv2LVasWFEhdSw/V1dXbNy4EX369IGDgwNyc3MrfJvFqUrnQWl4enpizpw5WL16dYnWz87ORs+ePXHr1i0AQFxcHHbu3ImuXbtWZDZJFadS2Rmoijp06ABNTU1cvXoV8+fPx6hRowAAYrEYJ06cwPLly3Hv3j04OTmBx+MBAHR1ddGzZ0/UqVOnRNsICgqCoaFhseuZm5tDR0cHAKCtrY2oqCioqamVcc8K3770diqTl5cXzp49i5MnT+Ljx48YP348AgIC0LFjR24dBwcH8Pl8TJ8+HdnZ2Rg8eDBq1qyJ8ePHc+s8ffoUbdu2xZAhQwAAenp6aNCgAff5p0+foKGhAU1NTWhoaCApKQkCgeDr7WglyX/c89fbOnXqICwsDDVq1KjQfFy8eBG+vr7YunUrhg0bBjs7O0yfPh2NGzeu0O0Wp6qcB6X1ww8/4N27d+jcuXOJ1ldTU0PPnj1lylsgECA5ObnceSnpta26+Vb3SwYjckVHRzM+n8/c3NwKfHbw4EHG5/PZkSNHypR2Tk4OmzFjRpm+u2TJEmZjY1Om70q8efOG7dmzp1xpVJTly5ezJUuWFPp5WFgYMzIyYh8/fuSW/fPPP6xTp04sKyuLW7ZkyRKWnZ1daDpr165l0dHR3N9OTk7MwsKinLmv+qZOnVrsOoqoY8WxsbFhO3furNBtfG/GjRvHli9fXubvu7m5MT6fX648pKWlsfnz55crjaooOjqarV27trKzUeGoy6AMJk+eDG1tbRw7dgwikahU3xWJRPj9998RHBxcpm0rKyuX6XsSnz59wqJFi8Cq6JxWnz59KnIfX758CbFYjLp163LL2rdvj9TUVLx58wYA4O/vjxs3bmDSpEk4f/48xGKxTBru7u44e/aszDIlpW//VNi3bx+ePn1a7HrlrWMlUdxxJqVX3jpc3uORnZ2NxYsXIykpqVzpVDWfP3+Gvb09cnJyKjsrFe7bvwpWADU1NXTp0gUJCQl4+/YtPn/+jCNHjsDS0hIfPnwAkNf8tn79euzatQs2NjZcU97du3fx9u1bpKSkwMHBAR4eHvDy8oK9vT3WrFmD48ePo0OHDrh8+TJevHgBe3t7uf2C9+/fh6WlJXr27IkTJ04AAN6/f4+5c+fC0NAQHz58gFgshqurK4yNjbFnzx6IRCK4uroiJSUFDx48gIODA6KiouRuRywW4+DBg9i6dSvmz5+PqVOnIjQ0FACQkpKCQ4cOwdLSEiEhIVi5ciVMTU0xc+bMIvuB09LSsHHjRmzduhWTJ0/G6tWrkZ6eDiDvB2L16tUIDQ3FmzdvsHr1aty9e7dAGp8/fwYAmYuOtrY2AHD94OHh4fj5558RERGB1atXw9bWFpmZmQDymv0k/aY7d+7Erl27ZNKPiorCb7/9BlNTUxw8eFDufmRmZsLNzQ1DhgzBs2fPsGrVKpiammLcuHGIjIzk1svKysKePXuwfv16DBkyBCtXrsSXL1+QkpKCw4cPw9LSEm/evMHQoUMxaNAgiEQiCIVCODk5YcOGDZg9ezYWLFjAlREA/Pvvv9i0aRPmzJmDQYMG4cGDBwCA58+fY968eVi9ejWePHmCfv36oWvXrnj06BGAvK6Y58+fA8jrcnF2dpZbbwtz5swZbN68GePGjcOkSZO4uiBPbGwsHBwcsHXrVowfPx5bt27lLqYBAQFYvXo1kpKScP/+faxevRovX76U+b5YLMbdu3cxdepU7N27F7dv34aVlRW6du2KO3fuICEhAXZ2dmjXrh3s7e0hFAq57zo7O2P16tVwdHTE0KFD4eXlBSAvkBw0aBAMDQ2xa9cu5ObmYurUqRg7diwXnOc/D8LDw7Fx40YMHjwYHz9+hJ2dHUxNTbFixQqIxWIcOHAAnTt3Rp8+fRAYGAgAuH37Nnr06IEJEyYAAKKjo7Fo0SLunJQ+d8LCwvD777+jY8eOGD16NJKTk/HgwQP07dsXHTt2hLu7e5HH5N27d5g/fz42bdoEe3t7xMXFyXx+7do1LF26FFu2bMHgwYNx/fp1AEBubi6uXLmCMWPG4NKlS3LTjoiIwODBg2FoaIjNmzdzx8/Lywvt27fn0pLm7u6OyMhIhIeHw8HBgat70dHR2Lp1K5YuXYoBAwZw51VZyjcwMBBr1qyBnZ0dvL29YWVlhW7duuHAgQMyeSnsPCnsWpudnQ0HBwds374d8+bNw+TJk7nryZkzZxAfH483b97AwcEBvr6+cHBwgKGhIVe/rl27BjMzM6xYsaLI7RSVtyqhspsoqqqiugwYY2zbtm2Mz+ezmzdvssTERK4bQdIMfeDAAebs7MwYY0wkEsk01e7evVumefrdu3fs559/ZgMHDmTXr19njo6OzMvLiwUHB7NevXrJNAMuX76c9erVi+3du5c9e/aMzZs3j/H5fHbv3j3GGGOenp4y+WCMsV69erHdu3dzf1tYWMj8LW8727dvZ1u2bOH+3rJlC+vWrRv7/PkzS0xMZCdPnmR8Pp9t3LiRhYSEsMePHzM+n8+uX78ut7xEIhEbM2YMe/DgAWMsr9tk9OjRBZqwbWxsimz2vHXrFuPz+ezSpUvcsoCAAMbn89n9+/dl1s3JyWH79u1jfD6fbdiwgVv+/PnzAmW0e/du1rFjR/bPP/+wT58+sR07djBjY2OWkpJSIA+pqanMw8OD8fl8tmjRIubp6cnu3LnDTE1N2fDhw5lYLGaMMfb777+z2NhY7jtdunRhq1atYvHx8Wz//v2Mz+ezTZs2sRs3bjAHBwfGGGNr1qzhunMEAgFr164dl/fw8HC2bds2Lh/79u1jbdq0YZGRkSwqKopZWFiwIUOGsIsXL7L4+Hg2ZcoU1r9/f279CxcuyDQJy6u3jOXVMekug4sXL7I7d+4wxv5Xl62srLj9lCYQCFjfvn1ZUFAQY4yxL1++MAsLiwLNrfnroDShUMh8fHxY69at2fTp09m9e/dYUlISmzhxIuvZsyc7cOAAi46O5uq6pM55eXkxPp/PdRUtWrSIjRgxgks3OjqamZqaMgcHB67bLi0tjfs8/3kQHR3NZs6cyTp16sTOnDnDYmNjmaurK+Pz+Wzt2rXM09OTJSQksKFDh7Lp06dz6SxdulSm/KTPyYSEBPb3338zPp/Pdu3axUJCQpifnx9r06YNs7OzY5cuXWJJSUlszZo1rEOHDiw3N1duGX348IF1796dxcTEcH+3bt2ay3tUVBQzMjJiYWFhjDHGduzYwbp06cIYYywzM5M9e/aM8fl8duHCBS7N/PUjODiYGRkZsZs3b3LLEhIS2MKFC+XmibGCdScnJ4ctXryY248XL14wPp/P3N3dy1S+4eHhbPTo0axfv37swIED7O3bt2z+/Pky539R50lh19o9e/awMWPGMMby6njPnj3Z1q1buTTyX5eioqIYn89nz58/55ZJd9kUtp2i8lYV0KDCMpIMJmSMQUdHB23btpX5PCEhAV5eXrCwsEDjxo1hZ2dXaFqtWrVCo0aNULduXVhbW8Pa2pr77McffyywfuPGjTFnzhwAQKdOnWBtbY2zZ8/C0tJSbrNhcU2JLVu2lNlOSkoKjh8/DhcXF27ZtGnTcOLECZw4cQJz585F8+bNAQATJ05E48aN0aJFC9SrV0/mDlnaw4cPERAQgJ9//hkAoKqqimnTpmHu3Lnw9vZGp06disyjRO/evdGwYUMcOnQIP//8MzQ0NLg7qSZNmsisq6qqitmzZyMrK4sbDKqqqlpo2urq6pgyZQoAYPDgwTh48CCioqJkuicAQEtLC6amptx6kpHZEydOxN9//413796hbt26ePDgAfT19bnvmZubIycnB7q6umjfvj2AvCchjIyM8OuvvyI6Ohrnzp3Dw4cPufxs2bKFGwx56NAhMMZw6NAhAEB6ejrMzMwQHR2N7t2744cffkDDhg0xbNgwAMAvv/yCzZs3F7q/8uqtPPv27cPw4cMRFhYGAGjQoAFEIhGSk5MLDMC7cOECeDwe+Hw+AEBDQwMTJ07E1q1bSzxoUFlZGSYmJtDW1kbbtm1haWkJAOjXrx8cHR0xY8YMAECjRo2gq6uLiIgIAHkDRydMmMANuq1fvz7+++8/Lt1GjRph8eLF2LhxIwQCAezs7GQGAec/Dxo1agQ+n4+goCCMGTMGADB8+HA4Ojqiffv23HHv0aMH1+oEFDzfpP+uX78+TExMuLQaNWoEADAyMkK9evUwdOhQAICVlRXc3NyQlJQkU4ck9u7diw4dOnD5bdiwIYyNjbnPa9eujaFDh3LlXb9+faSkpAAAatasWaLBhy1btkSvXr1w9uxZ9OvXDwBw9epVrn6VxLVr1xAXF4d//vkHQF7rT9euXZGYmFim8m3atCmaNm2KsLAwrh78/vvvePz4Mc6ePQsLC4tizxN519q0tDTUr1+fy7d0eckjuf5Lkz7OhV3TV61aVWje8l+/KgMFBGX06dMnAOBOuPwXgd9++w03b96EtbU1RowYAXt7+yLTU1JSkvuEQnH9eioqKujWrRtevHhRmuwXuZ3Xr18jNzdXJj86Ojpo0KAB/Pz8uPwCsieGurp6oV0GXl5eUFdXlymn1q1bAwD8/PxKHBCoqanB2dkZmzZtwtSpU9GxY0fEx8eDz+ejRYsWcr8zffp0HD58GMnJyXIvrhLSeZOMss/Ozpa7rmS/a9asyS3r0qUL/v77b0RFRSExMRFqamqFBoKSbdWuXZtb5u/vD8aYTADy66+/cv8ODg7GlClTMGDAgELTlN6HWrVqFfsoX3HBYmZmJqKjozFq1Cjo6uoWuS6Qd5zz1+PWrVtDLBbj7du3pXqKQEVF9vIkXdYSampqXJdB06ZNsWrVKly7dg3BwcGIiIgoMFZm3LhxuHLlCgICArigTlr+8620eSgJeed0/nQlQU1hx+/Ro0fcEzQS0k+G1KtXD5s3b8ajR4/w6tUrxMTEyJRFSccLTJo0CZMnT8b79+/RsmVLPH/+HBMnTizRd4G8OtusWTOZ82DmzJncv8tSvjweT2Y9TU1NGBsbIzo6mttmcedJ/jrap08fCAQCHD9+HOnp6cjMzCz3OCt52ykub5WNxhCUQW5uLry8vKCvrw8jIyO567Ro0QLXr1/H+PHjceHCBQwZMkTh7yiQ0NDQQK1atRSWnuRESExMlFmuq6tb4AQuTZqpqakyJ7YkIi9tmo0bN8bff/8Nd3d3LFq0CK9evSryufk6depAX18f9erVK/E2JD/4+QckFkVyt6ympoacnBzExMQgNTVVZp2iHuuS/AhI7nglJO9hyMnJgb+/f4HvKeJRscJI+o7fvn0rszwjIwNZWVkF1meMFRhUVtbjXFrp6emYNGkSVFVVsXDhQrRp06bAOjk5OdDT00NkZCROnz5dofmpSBkZGdx4Gnlyc3Nhb2+PyMhILFq0CN27dy/Tdrp27Qo+nw9nZ2cEBgaCz+eXavBiTk5OgboDKL7Oamtrc+dPWc6TwMBATJkyBT179sS8efNkWgsUqTLO4dKggKAMXF1dkZCQgNmzZxcaad++fRuamppYuXIl3Nzc8OXLF9y4cQOA/Oam8oiMjOROeEmTuPSdrVgsLtUPW+vWraGsrAxvb2+Z5SkpKejSpUuZ8mhiYgKxWCwzgEzyY1nWNIG85ux27dpxTa3yxMTEoEePHtwFQ5HlL12uiYmJUFFRQbt27dCiRQvk5OQUGOx0/vz5QtNq2bIlAMj8UDHGuC6Rli1bws3NTWbw2OvXr7mm/OKUZb+1tLSgq6uL/fv3y9ytnjt3Tm56JiYm+PjxI3e3BuQdZ1VVVZibm5d6+6Vx5MgRZGRkwMrKqtB19u7dC3t7e8yYMQM7duwoMBBPEVRVVWWCJUkdKe8dp7RmzZrB29u7QMuEZFvu7u7w9vYu1d18YSZOnIjLly/j2LFjRZ5nQME61rJlS/j7++P+/fvcsrS0NNy+fbtcecr/dFdiYiI6dOjAbbO054mjoyPMzMwKbWXMT3KdzX+ci7vOlvccrmgUEBRC3t0PkHdB3759OyZMmICxY8dyyyUVQfJ/T09P7sfP2NgYP/30E5o2bQogryk3OTkZSUlJ3EhcsVgs97EWkUhUoPInJSXJjJoPCQnB1KlTAeTdPfN4PFy4cAFBQUH4+++/kZmZiZCQEMTExADIa9qPiIjAx48f8e7duwLb+eGHHzBixAicO3eO+9EODAyESCTCiBEjAIC7EOW/yBX2GOavv/6KFi1a4OjRo1wZSbpUpF/2kZWVxe1bcc6dO4eoqCjs2LGDWyYWi7FgwQJuFHRSUhL27NmDJUuWcOtIWlPCwsJw//59iMViiEQiuSdzcSd4UFAQ9+8rV65g1KhRqF+/Ppo3b46+ffvi2LFjWLx4MU6fPo25c+dyd62SdKWPeZMmTdCvXz+cPXsWf/zxBy5fvgx7e3uufGxtbZGRkYExY8bgyJEjOHjwII4dOwYzMzMAeXeF8vIrOSaS/Q4NDeUu0PnrrWR96eNoZ2cHHx8f2NjYwNXVFZs2bUJ6errclxf99ttv0NbWxuHDh7llN2/exKRJk2S6QkpynPPnQ15eGWPcOjk5OQgNDcXr16/x77//4tmzZ0hPT8fr16+RmpqKp0+fol69ejAwMMCMGTOgo6MDBweHIrdZ2jwAeedgUFAQvL298ezZM1y4cAEA8ObNG3z+/FluGvl/TCTnVWH1b9KkSYiOjsamTZuQkZGBkJAQREVFISIiApGRkcjOzkZqairu3bsHHx8f3Lx5E0Dej090dDSX3/z7mn8ZkDdORnLNkIwdKkytWrUQExODz58/4+nTpxg0aBAaNGiAhQsXYvv27Th58iTmzp2Lvn37lrl8AXD7COQ9xRAcHIzJkycDKP48kXetzc7OxrNnzxAeHg53d3dER0cjMTERnp6eAPKumZGRkUhJScGLFy9Qv359qKur48qVKwgJCcE///yD2NhYrvwL205xeatsFBDI8fr1a+zbtw8AcODAASxevBhr1qzBlClT8PjxYxw+fBhr1qzh1v/06RPc3NwAACdOnEBKSgrEYjGmTZsGBwcHbN68Gf3794eFhQUAwNraGvr6+pgwYQIaNWqEK1eu4N27d3jy5AkuXrzIpXvr1i0EBATgxYsX3CNjc+fOhYmJCSZNmoSVK1fi1KlTcHFx4S62urq6mDt3Lk6dOoX169dj0KBB0NPTg66uLncCTZo0CY8fP8aePXtgaGgodztr165Fv379MGXKFKxbtw4uLi44ceIE1NXV8fHjR25/nZ2dERsbi7NnzyIuLg4PHz6Ej49PgTJVUVHB0aNHoaKigkmTJsHR0RGfPn3C1q1bAeQ1gZ44cQIBAQF4/vw53N3dER8fXyCdpKQkXL58GXv37oWqqiqcnJxk+hOVlJSgrKyMlStXon///ti/fz9Wr14t82PUqlUr9O7dG8uXL4dQKMTr169x//59JCQk4NSpU4iNjcXJkycB5AWAkvEihdWVXbt2wcHBAbVq1cKqVau4zzZt2oTBgwfj/v37OHLkCPr06YNu3bohKioKZ86cAQDs2rVLpovgjz/+wIABA3Du3Dns27cPAwYM4N7SKHl8tGbNmti3bx98fHzg6OgIHo+Hmzdvcj9Cjx8/RnBwMK5duwYgbzBiRkYGunfvDhMTE9ja2kJTU1NuvfXy8oKXlxcCAwO5x6QmTJgAe3t7REdHY8+ePVBWVpbpB5amqakJZ2dnhIeHw87ODg4ODtDS0sLixYsB5DWN7t27F0lJSbhz5w6uX78u81ilhKurK+Lj4/Hw4UO8evVKJj8HDx5EcnIyjh8/jri4ONy/fx++vr4YN24cGjVqhDlz5uD169eYNm0axGIx/Pz8EBQUhMWLF3PN3WlpaahXrx4ePnyI1atXIzk5ucB54Ovri3v37iEhIQGnT59GcnIy9u/fDyBvsNy7d+/g5eWFu3fvIj4+Hq6urgCAMWPGwNjYGDNmzMCzZ88wYsQItGjRAunp6fjy5Qv3DgwXFxfEx8dz5//Lly/x8OFDRERE4NSpUwDyzi95zclDhw7F4sWLcePGDVhaWuLKlSto3rw5mjVrBoFAgEGDBqF9+/ZYsWIFrly5glmzZkFDQwN3796Fnp4e13J1/fp1hIaGIiwsDFevXgWQd81LS0vjtlWjRg2MGjWqRIMJR44cidzcXEyfPh1GRkZQV1fHkSNH0LZtW7i4uODixYtYtmwZdHR0yly+QF436datW7Fz505s3rwZBw4c4ManFHWeFHatnT9/PhISEjBjxgzUrl0b/fv3R1BQENfNNW7cOISHh2PZsmVcC+qqVavw4MEDLFy4EF27dkXTpk3RrFkzpKenF7qdovJWFfCYItuxCPlOfPjwAb/88gtOnDhR4tfFElJdOTg4YOnSpSV+NXtFWrFiBWJiYmSegiKKQS0EhBBCCvX582cwxqpEMEAqFj12SEgZSPo0S/O4GSHVyZYtW8AYQ0REBBYtWlTZ2eEIhcIqMTPmt4haCAgppaioKBw5cgRA3hMBVWWEMCGK5Ofnh+vXr2PEiBFVZpa/q1evwtvbG4GBgTJ980QxaAwBIYQQQqiFgBBCCCHfaEAQGxsrMzXx0KFDcfz48XKl+eTJE0yePBmGhobcY3uV6eXLl7Czs8PcuXNhZWUFQ0NDmXeZV3XJycn4+++/0bt372Jn2iuNGzdu4PXr1wpLj1Q9L1++5N4z8eXLF6xevRqdO3dG7969sXfv3gLPrH/+/BlLlizBH3/8gYULFyI8PLxAmi4uLli2bBlWrFiBc+fOFbptX19fbNq0CS4uLvD19S0ynwkJCZg3bx7MzMzQt29f7lFCaTExMViwYAE2bdqEJUuWICEhocA60dHRWL58OfdIXn7Xr1/HokWLsHbtWpl1QkNDcfbs2RK9EOnWrVvo2bNnoa/qLi2xWIzbt29j0qRJ2Lt3b7nTk5410tLSEv/++y8yMzPx999/o1WrVjA0NMT+/fshEAgA5P0G/P7779zjjtXJzZs3uZkm586dy82oCOS9L+Hy5cuwsLDA2LFj8eTJEwB5Y5oOHDiAGTNmYPr06TA1NYWhoaHcmXKL9NWmUfpKYmJi2Nq1a5lIJOKWbdu2jT158qTcaUtm9Ms/q97XFhISwiwsLFh8fDxjjDGxWMyOHTvGjI2NFbqdjx8/ss+fPys0TYmkpCR29OjRAjPtKcKePXsUcrxJHpFIxIKDgyt8O4GBgcWu8/DhQ3bs2DHu7xUrVrDNmzezK1euMHt7e8bn89nevXu5z8ViMRs7diw7c+YMY4wxf39/1rt3b5aens6tc+bMGZmZ7oYOHcrN7ijt2LFjbNasWUwgEJRof6ZNm8Z2797NPDw82KRJkwrM0pmRkcH69OnD1dVbt26x4cOHM6FQyK0THh7ODh48yAwNDeXODvno0SPWp08fbobHGTNmsOPHj3OfBwQEcDNpFuXt27ds5cqVcmewLA3JMRQKhSwgIICZmJgUOqtlaWVnZzNzc/MCM6RKZnxNTU2VWe7r68vWrVtXrm1+rbqfn4uLC+Pz+SwgIEDu51u2bOFmjmWMsV27drEVK1Zwxy8lJYXZ2NiwVatWlWq731QLgUgkwpIlSzBz5kyZ920vXboUPXr0KHf6klffyntD29fk4eGBNm3acJPN8Hg8TJ48WSH7KO3vv/+WeUGJImlra6NVq1YVkradnR2cnJwqbO6I782NGze4Sa0qyocPH2RePCNPXFwc9u/fz72O98uXLzA0NMSKFSswcOBA7Nq1C127dpW5w7958yb8/Py4iYBat26NGjVqwNnZGUDeGxN37tzJvYFTSUkJ1tbW2Lp1q8ydtZubG86fP48///wT6urqxe7P+/fvMWLECMybNw+DBw/GkSNH0KxZM5m8ubq6Ijs7mztvf/nlF4SFheHKlSvcOk2bNoWdnV2BWSUltm3bhv79+3PXpkGDBmH37t3IyMgAkDeLYt26deW2TkgzNjbGpk2byvWCnPT0dPz9998A8iZPkszgqChqamro06cPXrx4we0fAO51ypKXqkk8e/as3JMIfY26L4+kjhU2R426urrMZ25ubrCysuKOX926dbF9+/YST2Il8U0FBB4eHqhdu7bcKYO/JZLXbAYGBsosl7wJURHc3d25N6pVlNJMklIaampqsLa2VkhT5fcuKCgI69atq9BtfP78Gfb29nJf3S1t586d6Nu3L1dvxGKxzOvDgbxzQHrSn1u3bqFZs2Yyb7Ns3bo1N6+Il5cXUlJSZLrbWrdujaioKG4Smo8fP2LDhg3YsGFDiYIBIG8OCMmUwUDemzp//vnnAnmT3q6ysjIMDQ25vEmTdxMSGhqK9+/fc7OGSvL+5csXmS7NiRMnYvfu3dwkWRUhOzsbixcvLjCxlaLP8V9//RXZ2dncFOEA0K5dOwAoMD/CmzdvyvVK4K9R9xUlJyeHCzAlGjRoUOpu5G8qIDh58qTMRDmMMTx48AC2trbcj4O/vz9WrlyJadOmwd/fH0OHDoWZmVmRk87Ik5GRgePHj6NVq1bYs2cPgLzX2A4dOpSbvz08PBwbN27E4MGDER0dDRsbG5iamnKvRZb4999/sWnTJsyZMweDBg3CgwcPitz20KFDkZWVhbFjx+Lo0aPcM7n5L44BAQHYvHkzFixYgAEDBnDvVH/58iUWLlyItWvX4saNG+jZsyd69erFPcYTFBTEzT++c+dO7Nq1C0BeC8zhw4exceNGjBw5EnPmzEFcXBwyMjLg6uqKAQMG4Pnz59i7dy+6dOmCAQMGIDw8HP/99x8GDx4MU1NTbh5waUlJSbC1tUWHDh1gb28vc9FMSUnBjh07sGrVKgwaNAibN29Gbm4uYmNj4eTkhO7du+P9+/fo06cPpk+fzn2vU6dOcHd3lzsj3OfPn3Hs2DH07dsXgYGBmDVrFkxNTTFjxgyZC1pZty1NKBTCyckJGzZswOzZs7FgwQKZV/UWduyfP3+OefPmYfXq1Xjy5An69euHrl27cnNfAMCrV6+wfv16bN26Fd27d8fWrVuRnJwMJycnGBoacsfz8ePHsLCwwIQJE7jvnjlzBtu2bcOaNWvQpk0buXXuy5cvOH/+PNLT03Hp0iU4ODggJSUFhw8fhqWlJd68eYOhQ4di0KBBEIlE8PT0xKJFi7Bjxw4MHTq0QN/ts2fPsGLFCqxfvx6jRo3ipuw+c+YM4uPj8ebNGzg4OMjtm09NTcXVq1dlzm9NTc0C0+WKRCLuBwLI6/OXfm01kDcrZWhoKHJycri7P+l1JHfjAQEBAPImTdLX18d///2HGTNmYPjw4dwYhsLo6ekVuNuWzptIJMLbt28L3EHr6Ohw84sUR1JOReVdsqx+/frw8PCQm05cXBycnJzQrVs3AJA5nz09PbF9+3Z06NABY8eOLbTF0N3dHZGRkQgPD4eDg4NMPQWAS5cuoUePHrC0tJTJW2HXlMJ069YNdevW5eZmAPLGd2lra+Phw4dcUPnp0yfo6+tzAUl0dDS2bt2KpUuXYsCAATh48CD3/bt372Ljxo3YuHEjOnTogJMnT8qt+5J9v3nzJjZv3gxbW1uMGjUKb968gVgsxr1797jfmi1btsDMzAy3bt3CoUOHYGlpiZCQEKxcuRKmpqaYOXOmQt+lMGLECDx58gQjRoyQmUAu/29CscrYxVHlxMXFMT6fzx4+fMgtE4lEzMfHh7Vt25brx/r48SP77bffmKWlJTt58iSLj49nK1asYB06dJDpu5Pn+fPnjM/ns2fPnnHLevXqJdNH5uTkxCwsLBhjjMXGxrJFixaxTp06scOHD7PY2Fjm5OTEjIyMuP7/8PBwtm3bNu77+/btY23atGGRkZFF5uXu3bvM3Nyc8fl8ZmVlxR49eiTzeUpKikz/kbu7OzM0NGQvXrxgQUFBzMLCgg0fPpw5Ozszf39/ZmNjw4yMjLg+K8m+Svfv79mzh717944xxlhWVhYbOHAgmzRpEktNTWUXLlxgfD6fOTo6Ml9fXxYREcG6devGfvvtN+bq6soSExPZ3r17WatWrVhcXJzMNlasWMH1Dbdp04bNnTuX2+ayZcu4/t7o6GhmbGzM9u3bx2JiYti6desYn89nhw8fZm5ubszJyYn7Xnp6OuPz+ez69esFyi4xMZEdOnSI8fl8tn79evbq1St27tw5ZmxszObNm1fubUtbs2YN27NnD2OMMYFAwNq1a8c2bNjAGCv62EdFRTELCws2ZMgQdvHiRRYfH8+mTJnC+vfvz61vZWXF5e/du3dsy5YtjLG8es/n89mFCxe4dZcuXcpsbGwYY4x9+vSJDR8+nPvs+PHjRY6LkU4rPj6e7d+/n/H5fLZp0yZ248YN5uDgwDIzM1m7du24/vAzZ84wIyMjLn/e3t5swIABXF/3smXLmLm5OTfWx8bGhi1fvrzQPFy5coUZGhoW238/depU9vjxY+7vdu3aydQnxvLOUT6fz+Li4tj69esZn8+XGVMQFRXF+Hw+O3DgAGOMsZ49ezIbGxsWHx/PxGIx27RpEzM0NGQvX74sMi/5DRgwgL1//54xljeGhs/ns82bN8uss3TpUta6desC37WwsCjQF3/48GHG5/Nlxl5Ijv3atWtl1p07dy6bPn263HzFxsay33//nfH5fMYYY6mpqezOnTuMz+ezJUuWMB8fH/bu3TvWpk0bdvjw4UL3b/ny5Vwdk873+PHj2Z07d1hMTAwbMGAAmzNnDvd5YdeUoqxatYq1a9eOZWRkMMbyxhC4u7szPp/P7t27xxjLq9OSepCTk8MWL17McnNzGWOMvXjxgvH5fObu7s6ys7NZ7969ubTv3LnDXFxcuL/zn0deXl7M2dlZJi+dOnViaWlp7M2bN6xt27Zs9OjR7MGDB2zlypXMy8uLnTx5kvH5fLZx40YWEhLCjUWTd22SkFxPIyIi5H6+e/du9vz5c+7v3NxctmbNGsbn8xmfz2cLFy5ksbGxRZajPN/Mmwols85J+tWBvOYqExMTmSj8hx9+QKNGjcDj8TB+/HgAQL9+/XDx4kUkJiZCX1+/VNvN3yQm/be+vj6aNGmCWrVqYdq0aQDATbgTHR0NXV1dHDp0CIwx7s45PT0dZmZmiI6ORpMmTQrd7i+//IKrV69i8+bNuHHjBqZPn45x48Zh7dq1UFJSgqurK1JTU7l0MzMz0aVLF8TExMDc3Bw//vgj6tevz/XH/v777+jfvz/Onz8vM3GTRE5ODk6fPg01NTVuZKuhoSFSUlJQp04drmmuX79+aNu2LQCgY8eOSEpKwrhx4wDkTeq0e/duREdHQ09Pj0t7zpw5aNSoEXr16oX09HTs3bsXycnJCAsLQ0BAgEz/Z/fu3ZGeno4ff/yRmzlwzJgxBV6rWrt2bdSuXRuBgYGwtraW+UxHR4fLo42NDVq0aIEOHTrAy8sLly9fRmpqKkJCQsq8bYno6GicO3eOa95UV1fHli1b0KBBAwAo8th3794dP/zwAxo2bMhNKvPLL79g8+bNXPoJCQnYs2cPFi9ejFatWnF3MPKaaaWXJSUlITAwEFevXsXAgQMxfPjwEj/poauri/bt2wMAhg0bBiMjI/z6668QCoWwtrbmmq/r168PsViMtLQ01K5dG7t375bp65a0ypS0STkoKAgaGhpFNtm/ffsWNWvWRM+ePbllPB6vQHO7ZCY9FRUV7i5euqVB8pSCiooKvnz5gri4OEyePJm7tsybNw+nT58u1Sx1d+7cQdeuXWFgYCCzPH/eRCIRN7VucYrLuzR9fX2u1S8/fX19GBkZcX9raWlxLyIaPnw4TExMAAB8Pp+bya80OnbsiD59+gAAfv75Z+58KOqaIhaLC60b1tbWOH/+PB49esSNv7C2tsaGDRtw69YtWFpawtvbm7vuXLt2DXFxcfjnn38A5B3/rl27IjExERkZGfj48SP++ecfTJ48GZaWlkU+pbRv3z60atWKO2dr1qyJVq1a4ePHj2jfvj20tbXRrl079O7dG7179wbwv2MyceJENG7cGC1atEC9evWKLMuSjOWQXkdFRQUbNmyAlZUVNmzYgGvXruHJkyf466+/ZM6H4nwzAYFkml55gzDyD6xQUlKSqWyS75S0Cac0A2+K21ZwcDCmTJkid/BLTk6OTJ6UlZVlTn59fX04OTlh9OjRWL16NU6dOgUjIyOMGTMGwcHBaNeuHezs7ArdB+m0mjVrhgYNGiAqKkru+lFRUUhPT8f06dPl7r+8wSv5m3MlPwZFve73559/xt69e/HhwwcEBQVBV1dXZh+k/y0p18J+kGvVqoXExES5n0n2QfqC3KVLF1y+fBlRUVHl3jaQ1z3FGJNp0v3111+5fxd17CXbyF93pOvDvHnzsHXrVty+fRuzZ8/mBsYVx9jYGBYWFli8eDFcXFywaNGiUk3QJMlT7dq1uWUqKirYvHkz3rx5g2PHjnHdIpIfXz8/P5nArGnTptx04CWRmppa6AArIO9cOXDgADZt2iSzXFdXt0DfuUAggLKyMurWrcv9yKenp3M3DpIBa9ra2ty/pQOR2rVro23bttwslZmZmTLT9daoUUPmB/nz58+4fPmyzDTd9erVg6qqaoG8ZWRklHggnnTepb8vybs0dXV1ubMmSuQ/fyXnh3T9U1dXL1Mzt3QaNWvW5Pq5i7umFKZLly6oV68ebt26hdzcXFhYWEBNTQ2Wlpa4f/8+4uLioKmpyQVWwcHBaNasmcz5Kz1b58SJE7F161ZcuHAB8+fPh5WVVaHbDg4OxsKFC7mgWN6+5r8mSPZfeh+LK0vJdwq7VgqFQrnjSnr27ImrV6/iwIED+PvvvzF//nzcvXu3QH0odLslWqsakPz4fI13XOf/oSuPnJwcbvCStOTkZBw8eBAdOnTg/pP0U1+6dElm3W7duuHkyZNQVVXl5rkvKt3CaGtrcz/a8vKZnZ2NkJAQmeUpKSkles65pDQ0NAD878cvKCiowDEtah+k8Xi8Uo2ylfS/qqmplXvbknQAyExvDID7ESjLMZI2ZcoUuLm54YcffsCaNWuwfPnyEn2Px+Nh7969+Ouvv5CQkICJEycW+ex9SW3cuBG3b9/G/PnzZQIfIK8s8pcDY4x7brw4NWrUKPLc3rlzJ+bOnQstLS2Z5UZGRgUGuiUmJqJ169ZQUlLi7oyly1zyLgBJ66KamlqBsSg6OjrctgYMGCBznkqmaQby7g63b9+ONWvWyJxbPB4PfD5fbt4kd+TFkeRdOg1JACxpAZPeXkUN4i0NHo/HXS/Kek1RUVFB37598fDhQ1y/fh2//PILgLxgOy0tDRs3buSWSbbz9u3bAulIjvnq1atx9OhRKCkpYd68eXBycip02/LSEovF3A2pokjqVkpKitzPk5KSZG40pH8T1NTUYG9vjwULFkAgEMDb27vE2638GqIgkqZ+eXOrl9eHDx9k7gCkm/JVVVWRlZXF/S0Wi2XWLU7Lli3h5uYmM5Dm9evXCAsLw4gRI+Dq6sr9J2nKj4yM5LpIJBo2bIhmzZqhfv36XLp3796VGcATHR0t82hO/he4JCYmokOHDgAKtoL89NNPUFVV5QZQSpw7d06hc3lHRERAX18fLVq0gIGBARISEnDmzBnu85ycHJkLblEEAgHXPF8Y6WOVkJAALS0tGBgYlHvbQN4xAPLmO5BgjMHd3Z37vLBjXxK3b9+GiYkJTp06hZkzZ3LdHUDR9fL9+/cIDw/HgAEDcO3aNXTt2hUnTpwo8X7J4+npCRcXF8ybN09uEGZgYIArV67I3BHfuHFDJo9F0dfXh0AgkPtDcezYMXTt2lXmffuSH8nBgwcjMDBQ5gmGoKAg7gmAHj16QFtbG//995/M5wYGBmjRogVUVVXRo0cPvHnzRmabnz9/RteuXQEAu3btkjlPe/Xqxa23Y8cOjBs3TqYrUjpvPj4+3HKhUIjQ0FCZpxOKwufzYWRkJJNGUFAQ6tWrV6DFpyTnQnmV9jpQnmuKtbU1MjIyIBKJoKmpCSDvWNauXRtPnz6VeQS7ZcuW8Pf3526WACAtLQ23b99GamoqvL290aNHD1y6dAlDhgwp8lxo2bIljh49KlOPr1+/LvMYpCK0atUKysrKcgeYCoVCvH//XuZ36M6dOwXW6969OwBwvwkl8c0EBK1atYK6ujo+fvxY4DORSCTz45ebmyv3R1vesi9fvsDKygqHDh1CZmYmevXqJfNMcOPGjfH48WMEBQXB3d0dL168QHJyMt6+fQuRSAShUFjktmxtbZGRkYExY8bgyJEjOHjwINc3+eOPP8Lc3Jz7T3LBY4xh5cqVMj8kvr6++PDhA2xsbAAA48ePR40aNbg3hR07dgzr16+XiZxDQkK4C+zLly8hEokwcuRIAP/r2ggLC8P9+/dRq1YtjB8/Hrdu3cL06dNx+vRprFq1imvelOyP9L7mD44k28pfHpL+69zcXDg7O2PNmjVQUlJCt27d0Lp1a2zevBnr16/HqVOnMGvWLO5kl6Qj73G1z58/QyAQyIw4lyc4OJj799WrVzF9+nSoqKiUa9sSTZo0Qb9+/XD27Fn88ccfuHz5Muzt7bnjWNSxl5SHvLojqcsHDx7kml9//fVX6OjocBfHxo0b49atWwgJCcGpU6cQHByMmJgYhISEQCAQcJMz1apVC71790azZs0K3Q91dXWEh4cjICAAHz9+lLvvkny4u7vLTDwTHByMwMBA2NnZITk5Gba2tvDw8MCePXvw5s0b7lxSV1dHZGQkUlJSuKcPpJmbmyM3Nxfx8fEyy8+ePcsFvY8fP8bDhw9x8OBB3Lt3D0DeY4gtW7bk+q19fX2Rk5PD9S+rqqpixowZ3AVVKBTi+vXrWLJkCbeN+fPnw9vbm3vDYVxcHPfUEJB3Ny59nkr2aefOnWCMITExEY8fP8aDBw+wbds2rlVo5MiRYIxxd5y3b9+GgYGB3IAgNzdXbvPxvHnzcO/ePe6YXLp0CfPmzSvQivnhw4ciWx4kdSr/LJ7561/+mwhptWrVQkxMDD5//oynT59y6xeWhoaGRpHXlKJ07twZ2traMt1Qkm6Dnj17yuz/oEGD0KBBAyxcuBDbt2/HyZMnMXfuXPTt2xe5ubk4cOAAGGNQUVGBlZWVTFdW/ro/ffp0xMTEYMyYMXB2doaTkxO8vb25R93FYnGBa4KkLPMHs0WVpb6+PsaPH4+jR4/KdOOmp6dj9erVGDZsmEzQ5O/vj7/++ksmzStXrsDU1BSmpqbFlien1MMQq7BFixaxP//8U2bZqVOnmJGRERsyZAjz9fVlnp6erEePHszU1JRdvXqVRUVFsSVLljA+n882bNjAkpKSZL4vFovZvHnzWPfu3dmWLVtYYmKizOc+Pj6sd+/erHv37uzy5ctsz549bNKkSezOnTvMx8eHDRw4kLVq1Yq5urqy2NhYtnnzZsbn89miRYu4Efx37txh/fr1Y+3bt2ezZs0qkIf8/vrrL8bn81nr1q3ZhAkT2MyZM9nEiRMLjHp+8eIFGzp0KDMxMWE2NjYsKiqK+8zGxoYNHTqUbdiwgf35559s9uzZMm/kys3NZXZ2dqxLly7s1q1bjLG8N4Vt2LCBmZubsx49erBDhw4xxvJGTEtGKS9atIhFRESwhw8fsp49ezIzMzPm7u7O4uLimIODA+Pz+Wz+/Pnsw4cPLDs7mx07dowNGzaM2dvbs2XLlsk8JcJY3oj4GTNmMBMTE2Ztbc2NHA4ICGA2NjaMz+ezdevWcU8uSDx//px17dqV5eTkyC1D6Sccdu/ezZYuXcr27Nkj86a2sm5bWlpaGlu0aBFr164ds7KyYjdu3JD5vLBjf+PGDda+fXvWq1cv9ujRIxYUFMSmTp3K+Hw+279/PxMIBKxNmzbM2tqabd++nS1fvpy9efOGS/fBgwesS5cuzNLSkj179oytWrWKzZ49m3l6erI3b94wPp/PJk2axP766y+2YsUK7qkXefbs2cM6dOjAduzYwSIjI7k3Atrb27Pw8HDGWN5I7hkzZjBTU1O2aNEiFhISwjp16sSWLFnCje4+deoU69mzJ+vUqRNbv349y8zM5Lbx8OFD1rlzZzZt2jS5TxIIhULWq1cvmScIrl27xgwNDbmR1ZL/WrVqJXOexsbGsgULFrCtW7eyJUuWyJwH0vu4du1atnz5cnb16tUCn3t6erKZM2eyvXv3spUrVxb79rojR44UyBefz2cdO3aUqZMhISFs3rx5bMuWLWzFihUsOTlZJp34+Hh27NgxZmRkxAYOHMhu375dYFunT59my5YtY2vWrJF5S6E0CwsLdvfuXbmfhYSEsAkTJnB169OnT8zR0VHmfL5x4wYzMzNjvXr1knnKStrbt29Zjx492KhRo1hCQgKX72HDhrEXL16w169fs2HDhjEjIyPuzZGFXVNKYsOGDQXepHr37l125cqVAusGBwez8ePHs7Zt27Jhw4YxPz8/xlhe+fL5fDZixAj2119/sWXLlrHQ0FDue9J1X+L06dOsd+/ezMzMjK1YsYIJBAImFArZiRMnmJGREevduzdX1jExMdz5smHDBvbp0yd25swZZmxszIYMGcL++++/QvdPJBKxo0ePsqFDh7IpU6YwW1tbNnHiRLlPJ/Ts2ZPx+XzWo0cPZmdnx6ZMmcLWrFnDUlJSSlyejDH2Tc12GBwcjBUrVtC0mCUwYcIENGzYEFu2bKnsrFSIbdu2QUdHB7a2tnI/9/LywsSJE3Hv3j00atToK+eOlIWrqysiIiJK/37271xISAgWLlyIy5cvK7R7j3x7vpkuAyCvT61Xr16FPl5Dvg9paWkICAjgHqkk34axY8ciKipKbrcgKdyJEyewceNGCgZIsb6pgADI61N7+/ZtoY/PkTxCofCrPJHxteXk5ODgwYPYunVrkc9z5+8rJVWfsrIyNm/eDGdn5xIPRvze3bx5E127di12LA0hwDcYECgpKWHRokUIDw8vctDG90okEuHMmTMIDAyEl5cXN/jqW+Hl5YVZs2bJvPgov8DAQO6FQ//88w9NglSNaGtrY+HChXj16lVlZ6XKi4yMxE8//VTgxVyEFOabGkNACCGEkLL55loICCGEEFJ6FBAQQgghhAICQgghhHxDkxuVh1AoxOfPn1GjRo0q8b5vQgghREIsFiM7OxtaWloFZrJUJAoIkPea2/wTrxBCCCFVSdOmTWVena9oFBDgf1PgNm3atMgpVstKJBIhODgYfD6/VLPvfeuoXOSjcpGPykU+Khf5vqVyyczMREREhNwpjxWJAgL8b+7pWrVqycx7riiS9yGoq6tX+4qpSFQu8lG5yEflIh+Vi3zfYrlUdJc2dZgTQgghhAICQgghhFBAQAghhBBQQEAIIYQQUEBACCGEEFBAQAipRujFYYRUHDq7CCFVXq5IDGVlZZiampb7EbJckVhBuSLk20LvISCEVHmqykoY5fwCCSmfoaFRGzwer0zpqKspw22iuYJzR8i3gQICQki1kJEjgiBHBJ6aqMwBASGkcNRlQAghhBAKCAghhBBCAQEhhBBCQAEBIYQQQkABASGEEEJQyU8ZHD16FIcOHULdunXh6OiILl26IDk5Gfv374eenh40NDQwfvx4AHnzQTs5OUFXVxe5ubmYNWsWAEAsFsPJyQlaWlpITU3F/PnzoaJCD08QQgghpVFpLQSBgYHQ1NTE06dPMWbMGCxcuBBisRhr167F2LFjYWdnB39/f4SEhAAAdu7ciW7dumHatGkQCoV49OgRAMDFxQX6+vqwtbVFs2bN4ObmVlm7RAghhFRblRYQ1K5dG6NGjYKqqiqmTJkCoVCI5ORk+Pj4wMDAAABgZmYGV1dXCIVCXLx4EZ07dwYAmJubw9XVFQDg5uaGLl26AAA6duzILSeEEEJIyVVaQNCoUSPu3yKRCNra2vD394e2tja3XFdXFyEhIYiIiABjDDVr1gQA6OnpISQkBJmZmQgJCUH9+vW59SMiIiASib7uzhBCCCHVXJXobH/w4AGmTp2K9PR0aGlpcctVVVWRlJSEtLQ0meUqKircch6PB01NTW65UChEamoqdHR0Sp0PkUhUIcGEJE0KVGRRuchH5VKQsrIyGGMAwP2/LCTf/ZbKluqLfN9SuXytfaj0gCA7OxuvXr3C8uXL8eTJE2RnZ3OfZWZmok6dOtDS0pK7vG7dumCMIScnBzVq1EBWVhYAoE6dOmXKS3BwcPl2phh+fn4Vmn51ReUiH5VLHiUlJZiamiIjQwAAEAgEZU6LqeVNjOTr6wux+Nua5Ijqi3xULiVX6QGBs7MzZs2aBR6PByMjIyQkJHCfxcbGwtjYGI0bN4ZIJEJ2djZq1KiBuLg4GBsbo0aNGmjevDni4uLQpEkTxMbGwsDAAGpqamXKC5/Ph7q6uqJ2jSMSieDn54e2bduWe6a2bwmVi3xULvKpq2tAkJMGDQ2Nsk9upJpXniYmJorMWqWi+iLft1QuGRkZFX7DClRyQHD69GlYWVlBS0sLGRkZiIyMRMuWLREWFobmzZvDz88PNjY2UFNTg7W1Nby9vdGzZ0/4+Phg5MiRAIARI0bA09MTTZo0kVleFsrKyhVacSo6/eqKykU+KhdZkiCAx+OVOSCQfO9bLFeqL/J9C+XytfJfaQHBmTNn8Mcff0BVVRVAXteBi4sLHB0dceDAAfz0008wMzNDmzZtAAALFy7Ezp07ERoaCg0NDVhZWQEAJk6ciD///BPHjx9HWloa5syZU1m7RAghhFRblRYQjB07FmPHjpX7maOjY4FlmpqaWLduXYHlampqWLVqlcLzRwghhHxP6NXFhBBCCKGAgBBCCCEUEBBCCCEEFBAQQgghBBQQEEIIIQQUEBBCCCEEFBAQQgghBBQQEEIIIQQUEBBCCCEEFBAQQgghBBQQEEIIIQQUEBBCCCEEFBAQQgghBBQQEEK+I6rKPOSKxApLT5FpEVLZKm36Y0JI1ZIrEkNVWTH3CIpMS5FUlJSgqqyE0SdeIiNHVK601NWU4TbRXEE5I6TyUUBACAGA7+qHMiNHhIzc8u0nId8aCggIIRz6oSTk+1X12vQIIYQQ8tVRQEAIIYQQCggIIYQQQgEBIYQQQkABASGEEEJAAQEhhBBCQAEBIYQQQkABASGEEEJAAQEhhBBCQAEBIYRUe0pKdCkn5Ue1iBBCqoiyzJ6orKwMU1NTKCsrlzst8n2juQwIIaSKKMsEU4wxCARfoKFRGzweD0D1mGCKVD0UEBBCFEpVmVdlpz+uDko7wRRjDIIcEXhqIi4gIKQsKCAghCiUipKSwqZSBgAdDVU4/9ZBATkjhBSFAgJCSIVQ1FTKtXKUi1+JEFJulRoQCAQCHDt2DCkpKVi7di233NbWFk+fPgUAXL58GYaGhsjMzISTkxN0dXWRm5uLWbNmAQDEYjGcnJygpaWF1NRUzJ8/HyoqFOcQQgghpVGpnXwCgQBCoRACgYBbFhQUhMGDB+Pp06fw9PSEoaEhAGDnzp3o1q0bpk2bBqFQiEePHgEAXFxcoK+vD1tbWzRr1gxubm6Vsi+EEEJIdVapAYGenh4aN24ss+z48ePw9/dHUlISdHR0AABCoRAXL15E586dAQDm5uZwdXUFALi5uaFLly4AgI4dO3LLCSGEEFJylT4MOP+o2GbNmiE6OhqjRo3CtWvXAAARERFgjKFmzZoA8gKJkJAQZGZmIiQkBPXr1wcA6OrqIiIiAiJR+fstCalM9KIZQsjXVuU62+3s7AAAd+/exerVq2FhYYG0tDRoaWlx66ioqCApKQlpaWng8XjQ1NTklguFQqSmpnKtC6UhEokqJJiQpEmBiqzvsVzE4BX7OJ7kRTMlkSsSQwlMEVmDsrIyGGNgrHzpSb6viLTypyf9/8rOm+T7iqy/ZTkG8sqlIvJW3XxL15evtQ9VLiCQ6NOnD9zd3REaGgotLS1kZ2dzn2VmZqJOnTqoW7cuGGPIyclBjRo1kJWVBQCoU6dOmbYZHByskLwXxs/Pr0LTr66+l3JRUlKCqakpBux/CEGOsNzpaaip4Nrs3njz5g3E4vK9lU6SN4HgCwTlfFSwFtQAAIIMAQTZ5d9PSXoZmRl56UqNOarMvDG1vKcffH19y13+QPmPgXS5KDpv1dn3cn1RhCobEABA48aNoa2tDV1dXYhEImRnZ6NGjRqIi4uDsbExatSogebNmyMuLg5NmjRBbGwsDAwMoKamVqbt8fl8qKurK3gv8qI7Pz8/tG3btsDrRb9n3225qNYED4Vf8PPePCeAhoZG0S+aUc0rMxMTE4VlTUOjNnhq5QsI1NXzzj8NdQ3wVMt/ZyNJT72WOgTZacWXy1fKm3oFlD9Q+mMgr75UVN6qk2/p+pKRkVHhN6xAFQgIpJu5MjIy8OHDB/D5fGRkZEBDQwMNGzYEAFhbW8Pb2xs9e/aEj48PRo4cCQAYMWIEPD090aRJE5nlZaGsrFyhFaei06+uvrdy4fF4JfpBK249yWeKLLuS5q24NBSVVv70ypuuIvNWEeUvSbcseZP+XkXlrTr6Fq4vXyv/lRoQxMfH48WLFwgODkZoaCgYY7C1tYWpqSnat2+PadOmcesuXLgQO3fuRGhoKDQ0NGBlZQUAmDhxIv78808cP34caWlpmDNnTmXtDiGEEFJtVWpAoKenhy1btsgse/Lkidx1NTU1sW7dugLL1dTUsGrVqgrJHyGEEPK9oGebCCGEEEIBASGEEEIoICCEEEIIKCAghBBCCCggIIQQQggoICCEEEIIKCAghBBCCCggIIQQQggoICCEEEIIKCAghBBCCMoREJw4caLAsvDwcFy8eLFcGSKEEELI11fquQyuX7+OnJwceHp6QlNTU+az1NRUHD58GMOHD1dYBgkhhBBS8UodEHTp0gXLli1DaGgoBAKBzGdqampYuHChwjJHCCGEkK+j1AGBtrY29u/fjwcPHqBfv34FPheJRArJGCGEEEK+njJNf6ympoZ+/fohKysLqampEIvFAADGGC5fvoxZs2YpNJOEEEIIqVhlCggAYMeOHfjnn38KtAjweDwKCAghhJBqpswBwZkzZ3DmzBkYGxtDWVmZW+7u7q6IfBFCCCHkKyrzY4cjR46Ejo6OTDAAAObm5uXOFCGEEEK+rjK3EGhoaGDVqlUFAgAvLy+4uLiUO2OEEEII+XrKHBAEBwdDXV0dHz58AI/HAwCIxWLEx8crLHOEEEJKT1WZh1yRGKrKinkZrSLTIlVXmQOC2bNno2XLljJdBowxBAUFKSRjhBBCykZFSQmqykoYfeIlMnLK9yi4upoy3CZSV/D3oMwBgaamJuLi4mSWJSYmwsvLC0ZGRuXOGCGEkPLJyBEhI7fqvRuGWi+qpjIHBJaWluDxeGCMcct4PB709PQwffp0hWSOEELIt4daL6qmMgcEjo6OGDNmjMyyq1evokWLFuXOFCGEkG9bVW29+J6VuZ0lfzAAAH379sWyZcvKlSFCSMlIBo4RQogilLmFYO/evTJ/i8Vi+Pn5FZjwiBBSMRQ5cExHQxXOv3VQUM4IIdVRmQOCK1euwNTUlHvkkMfjoVWrVtiwYYPCMkcIKZ4iml5r5SgXvxIh5JtW5oBg+/btMDExUWReCCGEEFJJyhwQmJiYwNfXF+fPn0dsbCyaNWuG0aNH06BCQgghpBoq86BCd3d32NjYICEhAU2bNoVYLIaDgwPu3bunyPwRQkiVRIM6ybemzC0Ely9fxp07d6Cvry+zfNu2bfjll1/KnTFCCKnKFDmoE6CBnaTylTkgaNeuXYFgAADS0tLKlSFCCKlOFPU8PQ3sJJWtzAFBdnY2/v33X7Rv3x5CoRBhYWG4cOECkpKSSpyGQCDAsWPHkJKSgrVr1wIAkpOTsX//fujp6UFDQwPjx48HAGRmZsLJyQm6urrIzc3FrFmzAOQ97ujk5AQtLS2kpqZi/vz5UFEp824RUgC9GpUQ8j0o1+RGDg4OsLW15V5fbGFhgT/++KPEaQgEAgiFQpl3F6xduxYLFy6EgYEBVq5ciZCQEBgYGGDnzp3o3r07evXqhT179uDRo0fo1asXXFxcoK+vj/Hjx+PixYtwc3PDuHHjyrpbhBRAz/oTQr4HpbrtcXV1xYkTJ+Dh4QFVVVX89ddfePr0Kc6cOQNXV1esXLkS2traJU5PT08PjRs35v5OTEyEj48PDAwMAABmZmZwdXWFUCjExYsX0blzZwCAubk5XF1dAQBubm7o0qULAKBjx47cckIUSdIsXK7/cmgAGiGk6ipVC0FOTg7Cw8OxZs0aqKmpAQC0tbW5IMDd3R0ikQhNmzYtcZqSFxsBwNu3b2UCCl1dXYSEhCAiIgKMMdSsWRNAXiAREhKCzMxMhISEoH79+tz6EREREIlEMtMyl5RIJIJIpPh3a0vSrIi0q7PqUi7KyspgjMlM5FUWku8Xl5b0eopIT5F5+9pp5U9P+v+VnbeK3M/SpCevXCpiPxV5nir6nJKXt+pyfSmJr7UPpQoI4uPj8fvvvxf6+dChQ7Fjxw4sXry4TJlJS0uDlpYW97eqqiqSkpIKLFdRUeGW83g8aGpqcsuFQiFSU1Oho6NT6u0HBweXKd8l5efnV6HpV1dVuVyUlJRgamoKgeALBOXsMqiFvCBakCGAIFtY7PrFvQa8tOlVx7Sk08vIzMhLtxyvR68O+1nW9KTLRZF5Y2p5N1e+vr4Qi8vfyqXIc6okeavK15eqplQBQUkG6yUnJ5c5M1paWsjOzub+zszMRJ06dQpdXrduXTDGkJOTgxo1aiArKwsAUKdOnTJtn8/nQ11dvcz5L4xIJIKfnx/atm1bppaLb1V1KhcNjdrgqZVzqlb1vIu0hroGeKqFp8UYg0AggIaGhkwLWlnTU2TevnZa0ump11KHIDut2HL5WnmrqP0sbXry6otC91M179xU9JtpFXJOFZG36nR9KU5GRkaF37ACZegyKE5kZGSZM2NkZISEhATu79jYWBgbG6Nx48YQiUTIzs5GjRo1EBcXB2NjY9SoUQPNmzdHXFwcmjRpgtjYWBgYGHDdGaWlrKxcoRWnotOvrqpDufB4vDL/CEmnUZq0iluvtOkpMm9fK6386ZU33eqyn2VJT/p7FbGfij5Hv1beqsP1pThfK/+lGlSYmZmJsLCwQj9/+vRpqQ+wdB+Snp4eWrZsyW3Dz88PI0aMgJqaGqytreHt7Q0A8PHxwciRIwEAI0aMgKenZ4HlhBBCCCm5UgUE06dPx8yZM3H48GGEhoYiKysLubm5iIyMxKFDh7BgwQLMmTOnxOnFx8fjxYsXCA4ORmhoKADA0dERzs7O+Oeff2BmZoY2bdoAABYuXIj79+/j+PHj0NDQgJWVFQBg4sSJCAsLw/HjxxETE4OJEyeWZpcIIYQQglJ2GTRu3Bi7du3CsmXL8Ndff8l8pqenh02bNnGPAJaEnp4etmzZIrPshx9+gKOjY4F1NTU1sW7dugLL1dTUsGrVqhJvkxBCCCEFlfrFRK1atcKVK1fw8uVLBAQEQCQS4aeffkL37t3L3HdPCCGEkMpV5jcVmpubw9zcXJF5IYQQQkgloRe0E0IIIYQCAkIIIYRQQEAIIYQQUEBACCGEEFBAQAghhBBQQEAIIYQQUEBACCGEEFBAQAghhBBQQEAIIYQQUEBACCGEEFBAQAghhBBQQEAIIYQQUEBACCGEEFBAQAghhBBQQEAIIYQQUEBACCGEEFBAQAghhBBQQEAIIYQQUEBACCGEEFBAQAghhBBQQEAIIYQQUEBACCGEEFBAQAghhBBQQEAIIYQQUEBACCGEEFBAQAghhBBQQEAIIYQQUEBACCGEEFBAQAghhBBQQEAIIYQQACqVnQF5GGPo378/wsLCoKamhkePHgEA9u/fDz09PWhoaGD8+PEAgMzMTDg5OUFXVxe5ubmYNWtWZWadEEIIqZaqZAvBw4cPsX79ejx9+hSPHz+GtrY21q5di7Fjx8LOzg7+/v4ICQkBAOzcuRPdunXDtGnTIBQKueCBEEIIISVXJQMCV1dX/Pvvv8jKykK9evWQmJgIHx8fGBgYAADMzMzg6uoKoVCIixcvonPnzgAAc3NzuLq6VmbWCSGEkGqpynUZZGdnw8TEBC9evICLiwv279+PrKwsaGtrc+vo6uoiJCQEERERYIyhZs2aAAA9PT2u5aAsRCIRRCJRufdBXrrS/yd5qku5KCsrgzEGxli50pF8v7i0pNdTRHqKzNvXTit/etL/r+y8VeR+liY9eeVSEfupyPNU0eeUvLxVl+tLSXytfahyAUGNGjVgb28PADh58iTWrFkDe3t7aGlpceuoqqoiKSkJaWlpMstVVFSQlJRU5m0HBweXPeMl4OfnV6HpV1cVUS48Hg88Hq/c6SgpKaFdu3YQCL5AkFO+k7IW1AAAggwBBNnCYtcXCAQKTa86piWdXkZmRl66xZTL18pbRe1nWdOTLhdF5o2pKQMAfH19IRaLy5UWkHdOmZqaKuScKkne6LpbclUuIJBmY2ODCxcuQCwWIzs7m1uemZmJOnXqQEtLS+7ysuLz+VBXVy9XnuURiUTw8/ND27ZtoaysrPD0q6uKLBcxeFBVVlyPmIZGbfDUynfxUlfPu0hrqGuAp1p4WowxCAQCaGhoFBnUlDQ9Rebta6clnZ56LXUIstOKLZevlbeK2s/Spievvih0P1Xzzk0TE5NypZOfQs6pIvL2LV13MzIyKvyGFajiAQEAGBgYwNzcHE5OTtyy2NhYGBsbo3HjxhCJRMjOzkaNGjUQFxcHY2PjMm9LWVm5QitORadfXVVEuSgDGH3iJTLKeQeio6EK5986KKTFQfL9kqZV3HqlTU+ReftaaeVPr7zpVpf9LEt60t+riP1U9Dn6tfL2LVx3v1b+q9ygwuTkZERFRQEAEhMTYWxsjIYNG6Jly5YICwsDkNcENGLECKipqcHa2hre3t4AAB8fH4wcObLS8k6qlowcETJyy/lfTvmbSAkhpDqoci0EQUFBWLx4MX7++We0a9cOEydOBAA4OjriwIED+Omnn2BmZoY2bdoAABYuXIidO3ciNDQUGhoasLKyqszsE0IIIdVSlQsIunbtCk9PzwLLf/jhBzg6OhZYrqmpiXXr1n2NrBFCCCHfrCrXZUAIIaTqUFXmIVdEXWffgyrXQkAIIaTqUFFSgqqykkIG6QL/G6hLqh4KCAghhBRLMki3vGrlVO8R/98y6jIghBBCCAUEhBBCCKGAgBBCCCGggIAQQgghoICAEEIIIaCAgBBCCCGggIAQQgghoICAEEIIIaCAgBBCyDdKSYl+4kqDSosQQki1VNQ8C8rKyjA1NYWycunejPg9z9tAry4mhBBSLRU1zwJjDALBF2ho1AaPxytReupqynCbaF4RWa0WKCAgVUKuSAxVZWqwIoSUnrx5FhhjEOSIwFMTlTgg+N5RQECqBJpNjRBCKhcFBKTKoNnUCCGk8lAbLalUNAqYEEKqBroakzIr72jcso4CJoQQonjUZUDKrLz9/pJRwE306uHEODMF544QQkhpUEBAyqU8/f6SUcAZOd/vc7+EEFJVUJcBIYQQQiggIIQQQggFBIQQQggBBQTfne/5Pd2EEEIKR4MKvzOKeiMgvQ2QEEK+LRQQfIcU8UZAehsgIYR8W6jLgBBCCCEUEFR11OdPCCHka6AugyqOZgEkhBDyNVBAUA3QLICEEEIqGnUZfCU0qx8hhJCq7JtoIYiIiMCpU6egqamJFi1awNraulLzkysSQ1X5fwGAZFY/QgghpKqq9gEBYwwrVqzAwYMHoaWlhenTp6NLly6oV69epeUpf7+/ZFY/DY3a4PF4JU6H+vwJIeTrUVXmFbihKw9FpvU1VPuAwM/PD2KxGFpaWgCANm3a4MKFC5g2bVql5ku6318yqx9PTVSqgID6/Akh5OtRUVJS2EBudTVluE00V1DOvo5vIiDQ0dHh/tbV1YWfn1+p0hCL8x7tEwgEEInKP3hPWVkZTesoISuXAQAYA7JU1VCzpjJKEQ9AqxYPGRkZMmmVhyLTU0RaknJpUPvb3s/SplXS+lLd97O06TWro4S6KP15VFF5qyrnp7z68i3uZ2nTKst1V5KeXg2GLKXy5a2mKkNGRoZCflOysrIA/O+3qqLwGGPlP8KV6MCBA4iIiMCWLVsAAOfOncPt27dx+PDhEqeRlJSEiIiICsohIYQQUn5NmzaVuQFWtGrfQqClpYXs7Gzu78zMTNSpU6fUaTRt2hQ1atSgpwEIIYRUKWKxGNnZ2VzXeEWp9gGBsbExrly5wv0dGxsLY2PjUqWhoqJSoVEXIYQQUh61a9eu8G1U+9vhdu3aISsrC1++fAEABAUFYdCgQZWcK0IIIaR6qfZjCAAgMDAQ586dg76+Ppo3b44+ffpUdpYIIYSQauWbCAgIIYQQUj7VvsuAEEIIIeVHAQEhhBBCKCAghBBCCAUEhBBCCAEFBIQQQggBBQSEEEIIAQUEFeLo0aPo3Lkz+vXrh+fPnwMAkpOTsXHjRhw6dAiurq6VnMPKIa9cAMDW1haGhoYwNDREUFBQJeaw8pw6dQrW1tbo168f3r59C4DqDCC/XACqMwBw5coVrFixAgDVFWnS5QJQXSkNCggULDAwEJqamnj69CnGjBmDhQsXQiwWY+3atRg7dizs7Ozg7++PkJCQys7qV1VYuQQFBWHw4MF4+vQpPD09YWhoWNlZ/eo+fPgAIyMj3LhxA7/88gv2798PAN99nSmsXKjOAOnp6Thy5Aj39/deVyTylwvVldKhgEDBateujVGjRkFVVRVTpkyBUChEcnIyfHx8YGBgAAAwMzP77qL4wsrl+PHj8Pf3R1JS0nc7n0SjRo3QoUMHAECHDh3QpEkTJCYmfvd1Rl65AKA6A+DEiRMYNmwYAFBdkSJdLgDVldKigEDBGjVqxP1bJBJBW1sb/v7+0NbW5pbr6up+dxG8vHLR0dFBs2bNEB0djVGjRuHatWuVmMPKJxaL8fLlS8yZMwdv37797uuMhHS5APju64y/vz9++uknaGpqAgDVlf+Xv1wAqiulVe1nO6zKHjx4gKlTpyI9PV1m2kpVVVUkJSVVYs4ql6RceDwe7OzsAAB3797F6tWrYWFhAXV19UrO4deXm5uL/fv348yZM0hPT0eXLl2ozqBgufzxxx/fdZ0Ri8W4efMmlixZgosXLwIA0tLSvvu6Iq9cAHzXdaUsqIWggmRnZ+PVq1cYPXo0tLS0kJ2dzX2WmZmJOnXqVGLuKo90uUjr06cPOnbsiNDQ0ErKWeVSVVXF/Pnz4ebmBg8PD9SqVYvqDAqWS3p6OvfZ91hnLl68iOHDh8sso+uL/HKR9j3WlbKgFoIK4uzsjFmzZoHH48HIyAgJCQncZ7GxsTA2Nq7E3FUe6XLJr3HjxjJNn98jPp8PAwMDtG3bluqMFEm5KCsryyz/3urMlStXuMGVAoEAOTk5SEhI+O7rirxyYYxh69at3DrfW10pCwoIKsDp06dhZWUFLS0tZGRkIDIyEi1btkRYWBiaN28OPz8/2NjYVHY2v7r85eLv74+6deuCz+cjIyMDGhoaaNiwYWVn86vLzc1Fbm4u1NXVIRAI0KpVK+jr63/3dUZeuQBAcHDwd1tnnJ2duX9fvHgR3t7e2LJlC+zs7L7ruiKvXBwcHL7rulIWFBAo2JkzZ/DHH39AVVUVQF4TuYuLCxwdHXHgwAH89NNPMDMzQ5s2bSo5p19XYeVia2sLU1NTtG/fHtOmTavkXFaOJ0+eYP369fj111+hra2N5cuXA8B3X2fklcvHjx+pzsjxvdcVeaiulB6PMcYqOxOEEEIIqVw0qJAQQgghFBAQQgghhAICQgghhIACAkIIIYSAAgJCCCGEgAICQgghhIACAkIIIYSAAgJSRQmFQpw9exYWFhaVnRWFOXz4MBwdHQHkzdNub2+Pffv2KXQbQqEQR44cwdKlSzF27FjEx8fLXS88PByLFy8u0fZnzJjBzRR3//59jBo1Cl5eXgrNtyJIly8hpPToTYWkQnh7e8Pd3R36+vpITU3F8uXLUbNmTe7zNWvW4Ny5cwAACwsLHDhwQOb7YrEYWlpa+Pjx41fNd0Xq1asX0tLSAABNmjTBly9fIBaLFbqN48eP44cffoCtrS1mzJiBgIAA6OnpFVivdu3aiI2NRdOmTYtNc9y4cWjZsiUAwNTUFGFhYQrNs6JIly8hpPQoICAKl5KSgpUrV+Lq1auoVasWnJ2dsWPHDqxevRoAkJCQgJo1a+LYsWMA8uYsl7h//z74fD4aNWqE1q1bV0r+FSUpKQleXl7o378/gLwJeiRq1aoFXV1dhW/z7t27mD9/Png8Hg4dOlToerq6uiV+r3uvXr24f9erV09mvvmqRLp8qwLpulxduLi4YMKECZWdDVJJqMuAKNzdu3dRt25d1KpVCwBgaWmJ06dPIyMjAwBw8uRJmJiYoGPHjujWrRt++OEHAHmztEk3+cqbEbG6yMnJwbJly5CVlVXoOhWxf4mJiVBSKtlpXdL18qvOx+VryV+Xq4Nnz57h5MmTlZ0NUokoICAK9+XLF8TFxXF/N2jQALm5uYiIiEBubi7u3LmDZcuWoVevXnj69Cm33pUrVxAbG4vDhw/j3r173PK3b99i0KBBsLS0xIcPHwps78aNGxg9ejRu3LiBcePGoVOnTtxUqEBef/nOnTuxcuVK2NraIikpCf7+/pg3bx527tyJqVOnYty4cQCA169fc+va29tzP+jy0ggMDMS8efOwd+9ebNmyBebm5ti+fTuAvLvD0NBQ3Lx5E8ePH0dkZCRWrlwJBweHQsvNw8MD+/btw/jx42XyLy0nJwe7du3Crl27YGdnh8OHDwMAMjMz8eeffyIlJQVubm5c60t+R48exY4dO7B+/Xr4+vpyy9+/f49Vq1bhwIEDGDt2LIKCggAAz58/x+TJk3Hp0qUCaYWGhqJv376wsrJCZGQkAODRo0cYPHgwYmNjZda9fv06hg8fjrt372LRokUwNzeHs7MzfHx8MHToUHTv3p3Lj0AgwKpVq3D48GFMnjwZV69e5dIwMzPDqFGjkJ2dDTs7Ozg6OiI7O7tA+Z46dQr9+vXDy5cvYWdnh44dO+LmzZt48OAB+vXrh759++LDhw9ISkrCmjVrYGlpCQAIDAyEtbU19uzZg7S0NOzcuRPDhw+Ht7c3hgwZAktLS4SFhWHfvn3o0aMHpk6dCqFQWKBs5NVlece3pPmMi4vDH3/8gZkzZ+Ls2bPo1KkTRowYgaioKAAAYwwuLi7Yu3cvRo8ejfPnzyMnJwdHjx5F7969cf/+fXTq1AmPHz/Gy5cv4eDggF27dmH8+PGIi4tDSkoKPDw8kJiYiD///BPe3t6YNWsW11rw8uVLdOvWDRcvXkRcXBw2btyImTNnYufOnejatSuSkpLkniOkmmGEKFhAQAAzMjJiV69eZYwx9vLlS8bn85mfnx+3TkJCAnNwcGDGxsYsICCAW87n81l0dDRjjLHo6GjG5/PZrVu3mFAoZLa2tszJyanA9jIyMlinTp2Yo6MjS01NZW5ubozP57MXL16w3NxcNnPmTJabm8sYY2zGjBls5cqVLDMzk02dOpX99ttv7P3798zDw4PFxsay3377jYnFYiYUClnnzp2Zh4dHoWlIlk+YMIFFR0czHx8f1qpVKyYQCBhjjNnY2LALFy4wxhjLzs5ma9euZcuXL+fyvXz5crZ7927GGGMvXrxg+/btY4wxlpiYyIyNjdmLFy8K7Ov27dvZ2bNnuf3u0qULu379Ove5hYUFe/78udzjcufOHbZkyRLGGGNisZgNHDiQ2/68efPYpUuXGGOMrV+/nv3++++MMcYyMzPZr7/+yu1H/m3cvXuXde/enQmFQsYYY0+fPmXXrl0rsO309HTWoUMH5uTkxDIzM9m1a9dYhw4d2K1bt5hYLGYbNmxgCxYsYIwx5uLiwlasWMEYY+z69ets4MCBXDoXLlxgHTt2ZLGxsczR0ZFbnr98JXXnzJkzLDc3lx08eJBZWFiwp0+fMsbyjuH27dsZY4z9+++/zMLCosBxyc3NZR4eHqxTp07My8uLicViZmtryyZMmMBCQ0OZQCBg3bt3Z48ePZJb3tJ1ubDjW9J8CoVCduTIEWZpacmePHnCEhMT2ZAhQ9ikSZMYY4y5u7tzx8/f358ZGxuzyMhI9vz5c8bn89nt27fZ1atX2YcPH9iwYcOYt7c3Y4yx6dOns6NHjzLGGHv+/LlMOZw/f57Z2Nhwf0vqc3Z2Ntu/fz/r0aMH8/HxYa6uriwnJ0fuOUKqFxpDQBTOyMgIO3bswMmTJ+Hl5QVVVVWoqKigSZMm3Dr169eHo6MjMjMz4erqig0bNhSanpWVFQCgdevWSEhIKPB5rVq1oKGhgX79+kFLSwujRo3CsWPH8PDhQygpKSEuLo5rCtXX1wePx0PNmjVRv359NGzYEAYGBjAwMMChQ4fQrl078Hg8KCsrw8PDA/Xq1YOvr6/cNFRUVKCpqQljY2M0atQIDRo0gEgkQmpqKtTV1WXyqKamhvr16xc6SNLDwwNKSko4fvw4AKBnz55ISUmRWUckEuHs2bNwdXXl9nvAgAFwc3ODtbV1UYcEAHDkyBHY2NgAyGv2lx6jMXPmTDRq1AhRUVGIiYmBjo4OAKBmzZrcv+WxtLTEtm3bcOvWLfTv3x8PHjzAsmXLCqxXu3ZtaGpqokuXLqhZsybatGmDL1++cMfWyMgIly9fBgD07dsX3bp1w+fPn/Hu3TsIBAIuneHDh+PixYuYMmWKTCtI/vKVjI/o3r07VFRU0LZtW5w5cwbdu3cHALRq1YprxSisC0RFRQV6enrQ0NBAp06dAABt2rRBbGwsmjdvDgD46aefkJycXGj5SBR2fM3MzEqUT2VlZdSrVw8//vgjevToAQCws7PDokWLkJWVBQ8PD7Rt2xbHjx+HWCxGly5dkJiYiB9//BEA0KdPH24/16xZgzZt2iAwMBApKSlcV15+hZWLmpoa9PT08NNPP8HExAQmJiZ4/fq13HOEVC8UEJAK0b9/f24w3fz58/Hzzz/LHYw2fvx47Nmzp0RpqqioQCQSlWhdAwMDpKWlITY2FvXq1cPkyZMLrMPj8WQuWjExMVBR+d8poa+vDwDFpiGdPwCFPjlQ1AXy06dPGDp0KAYOHAgAcreVnJyMtLQ0mSbqRo0a4dmzZ4WmKy0oKKhAoCLRoEEDHD58GO3atYOxsbFMk39R+ebxeJg4cSKcnZ3x888/Q1VVFWpqaoWuK5F//IKSkhLY/8/ErqOjAw8PD2hqaqJDhw7cI48S48ePx9KlS5GWlsYdo/zp589zUdsrSv50lJWVC6RTkidFSnJ8S5tPAwMDMMaQnp6OT58+Yd68eTA1NQUATJ06FQC4Ljbp/dDV1YWTkxN69uyJFi1alKgc8st/7hR1jpDqg8YQkAoVHByMx48fY/HixXI/z3+nWpySXrxycnLQrFkz6Orq4tWrV0hMTOQ+k+47l6anp4enT5/KbOPVq1elSqOsdHV1cfv2bZll+behra2NmjVrFnjsT/opjaLUrl0boaGhcj+bN28eunfvjj59+hT40SvOsGHDEBYWhk2bNpWopaI4+/btQ1ZWFsaOHVsggMnJyYGvry8mT56M9evXl+nHLD8lJaUSB5plVZLjW1q5ubmoXbs2dHV1C6SfnZ2NwMDAAt9hjGHSpEn47bff0LVr1yLTL025fI1zhFQ8CghIhUlLS8Pq1auxdetWGBgYAACio6O5C1dubi48PDxga2vLfUdVVRVpaWkICwuTe7Ev6gdA0p2QnZ2NsLAwDB48GO3atYOOjg7s7Ozw8OFDeHh4wM/Pj/uO9AXP2toa0dHRWLt2LXx9fbFz507UqVOnyDTEYnGBPEn+lt4XyXLpdaX/HjhwIG7duoWNGzfCy8sLmzdvRr169WTSVVZWxqhRo3DhwgVuma+vL3777Tfub6FQKHeQGwD8+uuvOHnyJGJjYyESifDp0yckJydDKBTi3bt3SE5OxufPn/H27VtkZWUhOjq62HwDgLq6OkaPHg0fHx+YmJjIPzhyvifvcwAICAhAcnIycnNz4e3tLZOXI0eOYPLkyZg9ezY+fPggUxbS6ef/f1Hb09HRQVJSEqKjo+Hv74/g4GCuXIrLc1HbkD7+hR3f0uQTgEyXmbe3N3fsBw4cCGdnZ+zduxdeXl7YtGmTzOOOknqempqKmJgYpKSkIC4uDiEhIVz5qqqq4suXLxAKhYiIiICOjg7Cw8ORmpqK58+fy9QX6TQBFHuekeqBAgKicAkJCfDw8MDevXvh6OjI9RMDec3ev//+O4YPH44tW7Zg2rRpMl0JgwcPxoIFC5Cens5d7M+dO4fo6Gi8fPkSfn5+hd7lPnjwAM7Ozti2bRscHR2ho6MDNTU17N+/H0pKSliyZAlevHiBUaNGwd/fH//99x8ePHjA3Uk1b94c27dvx7NnzzB37ly0bNkSfD6/0DRCQkLg4+MDb29vREdHw93dHQBw+fJl5ObmwtraGgcPHsTr168RExOD58+fc/kPCwvDf//9h3///ReRkZHo3r07Vq9ejZs3b2LFihUwMzND48aNC+zj4sWLoauri/nz5+Ovv/5C586d0bNnT4jFYpw7dw4JCQlwd3eXe3e4YMECmJubY8SIEVizZg1q164NgUCA6OhoTJkyBWvXrsXmzZthYWGBV69eISUlBb6+vnj//j0eP36MpKQkPH36FPHx8bhx4wY+f/7Mpd2vXz8MGzas0Drx6NEjJCQk4M6dO4iLi+PKys3NDbGxsbh37x5CQkLw+vVrjB07FtevX8eUKVPQvn175Obm4t9//8X58+dx48YNqKqqgjGGhg0bYsuWLbh7926B8j1//jyAvNH+cXFxuHnzJhISEnDv3j2Eh4fD09MTPj4+CAoKQvPmzTFw4ECMHDkSnp6eaN26Nb58+YKQkBDuew8ePEBkZCS3jYCAAHh5eSEkJAQPHjyQuTOWV5cLO76lySeQ10Kyd+9eHD16FB8/fsTcuXMBACNHjsS0adNw8uRJrFu3DoMGDYKamhr38q/Dhw8jNzcX9erVw/DhwzFt2jT8888/6N27N+7evQuxWAwjIyP8+OOPmDx5MrS0tNCtWzcYGRlh4MCB3LiJmJgYJCUl4c6dOwgMDMT9+/cBoNBzhFQvPKaINjdCKpmlpSU2b96Mzp07V3ZWvkv79u3D6NGjK+RlSyTPxYsXcenSJbi4uFR2Vsg3igYVkm9CSZp2iWIxxnDixAloa2sjKSmJgoEKRnWcVDTqMiDVnoeHBxISEnD9+nW5TbekYmRkZODo0aM4e/Ys13RNKkZUVBTu3LmD0NBQeHp6VnZ2yDeKugwIIYQQQi0EhBBCCKGAgBBCCCGggIAQQgghoICAEEIIIaCAgBBCCCGggIAQQgghoICAEEIIIaCAgBBCCCGggIAQQgghAP4P5xi+25qsk08AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 510x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams[\"font.family\"] = \"Times New Roman\"\n",
    "axs = plt.subplots(figsize=(5.1, 3))\n",
    "axs = sns.histplot(tmax95_series, bins=20)\n",
    "axs.set_xlabel(\"95th percentile of daily maximum temperature\", fontsize=10)\n",
    "axs.set_ylabel(\"Count\", fontsize=10)\n",
    "axs.set_title(\"Distribution of 95th percentile of maximum daily temperature \\n(in June-September) per census tract (2006-2010) in the Western US\", fontsize=12)\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"figures/95th_percentile_distribution.png\", dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b7197119",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GEOID\n",
       "04001942600    37.058506\n",
       "04001942700    34.900247\n",
       "04001944000    30.802767\n",
       "04001944100    33.009674\n",
       "04001944201    35.484156\n",
       "                 ...    \n",
       "56043000200    34.085806\n",
       "56043000301    36.470632\n",
       "56043000302    36.516460\n",
       "56045951100    35.442859\n",
       "56045951300    35.857422\n",
       "Length: 18108, dtype: float64"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmax95_series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "5a1484da",
   "metadata": {},
   "outputs": [],
   "source": [
    "tmax95_series = tmax95_series.reset_index(name='tmax95')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "5077cb0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove the census tract from Texas 48141010222\n",
    "\n",
    "tmax95_series = tmax95_series[tmax95_series['GEOID'] != '48141010222']\n",
    "\n",
    "# also 06013990000\n",
    "\n",
    "tmax95_series = tmax95_series[tmax95_series['GEOID'] != '06013990000']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "c2bdf1d3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GEOID     12943\n",
       "tmax95    12943\n",
       "dtype: int64"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmax95_series[tmax95_series['tmax95'] > 32.22].count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "793f367d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18106"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(tmax95_series)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86c0a12e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# In our study, 77% (n=13,579) of census tracts had 95th percentile summertime temperatures > 90°F (Fig. S5A). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "e3e770a0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GEOID     71.484591\n",
       "tmax95    71.484591\n",
       "dtype: float64"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmax95_series[tmax95_series['tmax95'] > 32.22].count()*100/len(tmax95_series)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "841acbbd-c6b3-4689-8d7a-8f4addd00055",
   "metadata": {},
   "source": [
    "## Wildfire day"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "efd3b9da",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "wfday\n",
       "0.0    99098100\n",
       "1.0      115632\n",
       "Name: count, dtype: int64"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['wfday'].value_counts() # western us"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "e53de808-278d-4e1c-b400-c9a205ed905c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df['wfday'] = np.where(df['wfday'] == 1.0, True, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "70353c82-ee62-4dfd-8284-7085e7bd5624",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "wfday\n",
       "False    99098100\n",
       "True       115632\n",
       "Name: count, dtype: int64"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['wfday'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "5ff7dbde-decf-4066-9597-b977ed72197a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = df.drop(columns=[\"tmax95\"]) #, \"tmin05\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c2ea9fe-3635-4b41-a88e-0bd96d1b97e7",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Smoke Polluted Day"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "61778296-c0f8-4a09-93de-36c0d81496fe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df['smoke_pm_non_zero'] = np.where(df['smoke_pm'] > 0.0, True, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "91a197e8-fac8-4ebc-a8a5-4e418b0887b8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# As seen in https://www.researchsquare.com/article/rs-2866201/v1\n",
    "\n",
    "df['smoke_pm_gt_five'] = np.where(df['smoke_pm'] > 5.0, True, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "701e5cf0",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['smoke_pm_gt_ten'] = np.where(df['smoke_pm'] > 10.0, True, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "27802353-675c-4699-a843-26b9cb701784",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time</th>\n",
       "      <th>GEOID</th>\n",
       "      <th>wfday</th>\n",
       "      <th>tmax</th>\n",
       "      <th>smoke_pm</th>\n",
       "      <th>heatday</th>\n",
       "      <th>smoke_pm_non_zero</th>\n",
       "      <th>smoke_pm_gt_five</th>\n",
       "      <th>smoke_pm_gt_ten</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2006-01-01</td>\n",
       "      <td>04001942600</td>\n",
       "      <td>False</td>\n",
       "      <td>15.978498</td>\n",
       "      <td>0.0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2006-01-01</td>\n",
       "      <td>04001942700</td>\n",
       "      <td>False</td>\n",
       "      <td>14.836937</td>\n",
       "      <td>0.0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2006-01-01</td>\n",
       "      <td>04001944000</td>\n",
       "      <td>False</td>\n",
       "      <td>14.192923</td>\n",
       "      <td>0.0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2006-01-01</td>\n",
       "      <td>04001944100</td>\n",
       "      <td>False</td>\n",
       "      <td>14.254541</td>\n",
       "      <td>0.0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2006-01-01</td>\n",
       "      <td>04001944201</td>\n",
       "      <td>False</td>\n",
       "      <td>15.200634</td>\n",
       "      <td>0.0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        time        GEOID  wfday       tmax  smoke_pm  heatday  \\\n",
       "0 2006-01-01  04001942600  False  15.978498       0.0    False   \n",
       "1 2006-01-01  04001942700  False  14.836937       0.0    False   \n",
       "2 2006-01-01  04001944000  False  14.192923       0.0    False   \n",
       "3 2006-01-01  04001944100  False  14.254541       0.0    False   \n",
       "4 2006-01-01  04001944201  False  15.200634       0.0    False   \n",
       "\n",
       "   smoke_pm_non_zero  smoke_pm_gt_five  smoke_pm_gt_ten  \n",
       "0              False             False            False  \n",
       "1              False             False            False  \n",
       "2              False             False            False  \n",
       "3              False             False            False  \n",
       "4              False             False            False  "
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "cce21b4b-98e0-4629-9041-335d6013599e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df.to_parquet(\"outputs/d1-heat-wf-smokeday-230929.parquet\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5d7955b-c15f-48cf-868f-1a00dd2a8fc3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
